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Abstract: Finite differences have been widely used in mathematical theory as well 
as in scientific and engineering computations. These concepts are constantly men- 
tioned in calculus. Most frequently-used difference formulas provide excellent approx- 
imations to various derivative functions, including those used in modeling important 
physical processes on uniform grids. However, our research reveals that difference 
approximations on uniform grids cannot be applied blindly on nonuniform grids, nor 
can difference formulas to form consistent approximations to second derivatives. At 
best, they may lose accuracy; at worst they are inconsistent. Detailed consistency 
and error analysis, together with simulated examples, will be given. 



1. Background. The fundamentals of calculus were developed over a very 
long period of time. They can be traced back to the late 1660s and 1670s, when 
Newton and Leibniz introduced the concept of derivatives through finite differences. 
The concept quickly became a cornerstone of calculus [2, 10, 13]. 




Figure 1.1. LEFT: Sir Isaac Newton [11, 13]. RIGHT: Gottfried Wilhelm von Leibniz [2, 12]. 
The two scientists founded calculus as well as the finite difference calculations. 

Different finite difference formulas have been constructed and utilized for ap- 
proximating the rates of change of a function in applications [2, 5, 6]. The function 
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is given either as an analytical expression or as a set of numbers at discrete points 
in the region of interest. The mathematical models of most science and engineering 
problems require such an approximation if computer simulations are desirable [3, 8]. 

In many discussions of one-dimensional finite differences, the aforementioned 
discrete points are uniformly distributed, that is, the distance between any two 
neighboring points is a positive constant h [3, 5, 6]. 

Let f(t) be a differentiable function on (a, b). Recall the definition of the deriva- 
tive at a fixed point t G (a, b) [10]: 

m = ,, m /('+*)-/('), 

J w h-*> h 
Thus, it is reasonable that the fraction 

9 (t,h) = f{t + h) h - f{t \ 0<h«l, (1.1) 

would provide an approximation of the derivative f'(t) on the set T = {t, t±h, t± 
2h, . . .} fl (a, b). In fact, it can be shown by Newton-Gregory interpolation formulas 
that g(t,h) is indeed an approximation of f'(t) [13]. Formula (1.1) is called a first 
order forward difference [3, 10]. Based on it, a second order forward difference could 
be constructed: 

gjt + h,h)- gjt, h) _ f{t + 2h) - 2f(t + h) + f{t) ^ q < /i -c 1. (i. 2 ) 
h h 2 ' 

We will leave the discussion of (1.2) to Section IV. 

This paper studies finite difference approximations on sets of discrete points, 
where distances between neighboring points vary. Many questions remain in the 
situation, such as: 

1. How good can the formula (1.1) be in applications? 

2. Are there different, or better, approximation formulas? 

3. Can the formulas be used repeatedly for approximating higher derivatives? 

4. How can we evaluate the errors in approximations? 

5. Can we demonstrate our study through the use of computer simulations? 

The preceding questions motivated our research in the subject. In this article, 
we will study some of the most popular finite differences for approximating the 
first and second order derivatives. Cases involving different types of grids over the 
intervals will be investigated. Section II will be for finite difference preliminaries. In 
Section III, we will concentrate on the first order differences. We will show that the 
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same difference formula may possess different accuracies on different meshes. Then, 
we will show in Section IV that, although second order differences are in general 
acceptable on uniform meshes, naive applications of the formulas on nonuniform 
meshes will lead to incorrect results. Computer simulated numerical examples will be 
given to illustrate our conclusions in Section V. Our investigation will be summarized 
in Section VI. 

2. Preliminaries. Let us consider an interval [a, b], where oo > b > a > — oo, 
defined in R. We will start with definitions of the uniform and nonuniform meshes, 
which are discrete sets of numbers used as domains for our finite differences. Meshes 
are also called grids in science and engineering computations [8]. 

Definition 2.1. A mesh over the interval [a, b] is defined as a set of m + 2 distinct 
numbers (m > 0) denoted as T = {to, t\, . . . , im+l}> in which tk+i > tk, k = 
0,1, ... ,m, and to = a, t m+ \ = b. The value tk is called the kth mesh point, or kth 
grid, of T. Moreover, we call the quantity 

hk = tk+i -t k , m>k>0, 

the kth step size of the mesh. Step sizes are often less than one. T is called a 
uniform mesh if hk = h > 0, k = 0, 1, . . . , m, otherwise T is a nonuniform mesh. 




Figure 2.1. Plot of the trigonometric functions with different frequencies on the uniform mesh 
(LEFT) and the nonuniform mesh (RIGHT). Blue squares or diamonds on the i-axis indicate 
the mesh points used in each case. While dotted curves are the same between the left and 
right figures, locations of the plotted function values (large dots on the curves) corresponding 
to different meshes are significantly different. 

Example 2.1. Consider meshes with twelve points over the interval [0, 1]. A uniform 
mesh is Tj = {t^ = k x /i}^.L with h = 1/11 (blue diamonds in Figure 2.1), 
while a nonuniform mesh can be T2 = {to = 0, tk+i = tk + hk, k = 0, 1, ... , 10}, 
where hk, k = 0, 1, . . . , 10, are calculated based on an equi-arclength formula for 
y = sin(27ri); that is, hk, A; = 1,2, ...,11, are obtained by evenly dividing the arc 
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of the curve of the function over [0, 1], and then vertically projecting the resulting 
arc pieces with equal arclength to the t-axis (the mesh steps are separated by blue 
squares in Figure 2.1). A composite trapezoidal rule is employed to evaluate the 
arclength. Let y\ = sin(27rt), yi = — 0.5 sin(47rt), and y-& = — 0.8sin(7rt) be three 
trigonometric functions with distinct frequencies. Their graphs over Ti (LEFT) and 
T2 (RIGHT) are given in Figure 2.1 as red, green and cyan dots, respectively. The 
true functions (continuous curves) are shown for comparisons. Matlab® is used 
[7]- 

Definition 2.2. Let function y = f(t) be defined on the interval [a, b] in addition to 
n > 1. We say that / is n times continuously differentiable on [a, b] if in additional to 
the continuity of the nth derivative f^ n \t), t £ (a, b), both directional derivatives, 
f( n \a + ), f^ n \b~), exist. We further say that / is sufficiently smooth on [a,b] if n 
can be as large as we wish. 

Definition 2.3. Let functions y = f{t) and y = g(t) be defined on the mesh T, 
where g is considered as an approximation of /. If 

\f(t)-g(t)\=0(hP), t€T, 
where h = max hf-, then we say that the approximation at t is accurate to the 

0<k<m 

order p with respect to the step sizes. From the approximation point of view, an 
approximation g is consistent if and only if p > [1,4]. 

Definition 2.4. Let / = {/ , fi, . . . , f m+ i}, g = {go, 9i, ■ ■ -,9m+i} be two func- 
tions defined on the mesh T, where g is considered as an approximation of /. We 
define the scaled local difference between the two functions at tk as 

sld(f,g) k = ^—^, k = 0, l,...,m + l, 
where a = max I fu\ > 0. The scaled qlobal error indicator is defined as 

sgei(/,sO = n max |sld(/, 5 ) fc | . 

Note that the above definitions are different from standard definitions of local and 
global relative errors, where signs are rarely considered. The value of sld(/, g)k offers 
not only scaled relative error information, but also the direction of the error, that 
is, whether g^ is greater or less than The latter is particularly useful if approxi- 
mations of oscillatory problems are investigated. The function sgei(/, g) provides a 
scaled overall error estimate and is easy to use for its simplicity in structure. Both 
definitions can be used when some of the f k values are zero. 
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3. First order differences. We assume that the function y = f(t) is n 
times continuously differentiable on [a, b]. Let T be a mesh over [a, b]. We define the 
forward, backward, and central differences as follows: 

D + f(t k ) = /(t ; +l) ~ /(tfc) , k = 0,l,...,m, (3.1) 

l k+l — l k 

D _ m = /fa)-/fa-i) > fc = lj2> ... >ro + lj (3. 2) 
tk - tk-1 

Sf(t k ) = Jfc = i >2 m. (3.3) 



Theorem 3.1. Let T be nonuniform and f be twice continuously differentiable. 
Then the forward, backward, and central differences are first order approximations 
of the derivative function f' onT. Further, 

D + f{t k )-f'(t k ) = (3.4) 
D-f{t k )-f'(t k ) = -^I/"( Cl ), (3.5) 

sf(t k )-f\t k ) = 2{hk l hk i) [hint*) - hurm , (3.6) 

w/iere > & > ifc, tk > Q > ife-i, ^ = 1, 2. 

Proof. Since the proofs of the forward and backward difference approximations are 
similar, we only need to show the cases involving forward and central differences. 
First, according to the Maclaurin series expansion [10], we have 

f(t k+ i) = f(t k ) + h k f>(t k ) + ^f"(t k ) + ... + -tf^(0, (3-7) 

h 2 h n 

/te-i) = f(tk)-h k - 1 f>(t k ) + ^f(t k )-... + (-ir^f^(o. (3.8) 

Therefore, letting n = 2, 

f(t k+1 )-f(t k ) = hk f(t k ) + %f"(0, 

f{t k+l ) - /(t fc _!) = (/l* + fcfc_i)/'(t fc ) + ^ - /ltl/"(0] • 

Equations (3.4) and (3.6) become obvious. Since h k /{h k + h k -\) < 1 and h k -±/ \h k + 
h k -i) < 1, according to Definition 2.3, the differences involved are indeed first order 
approximations. □ 
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Theorem 3.2. Let T be uniform. Then the central difference becomes a second 
order approximation of the derivative function f if f is three times continuously 
differentiate while the forward and backward differences remain as first order ap- 
proximations of f if f is at least twice continuously differentiable. Further, 

D + f(t k )-f'(t k ) = (3.9) 
D_f(t k )-f'(t k ) = ~f"(Ci), (3-10) 
Sf(t k )-nt k ) = ^ [/"'(&) + /"'(&)], (3.11) 
where t k+1 > & > t k , t k > Q > ifc-i, £=1,2. 

Proof. The corollary is a direct extension of the Theorem 3.1 due to h k = %_i = h, 
Definition 2.3, and (3.7), (3.8). □ 

Remark 3.1. Though D+f(t k ) = D-f(t k +i) for k = 0,1,..., m, on any mesh 
considered, the central difference is not a simple combination of the forward and 
backward differences on a nonuniform mesh. The central difference can be viewed 
as an arithmetic average of the forward and backward differences only on any uni- 
form mesh. Interestingly, this average improves the accuracy of the approximation 
significantly. 

Remark 3.2. The right-hand-side of (3.4)-(3.6) and (3.9)-(3.11) can be viewed as 
errors of the respective approximations. 



4. Second order differences. Let function y = f(t) be n times contin- 
uously differentiable on [a, b), where n is sufficiently large. We will investigate if 
we can approximate the second order derivative function by applying the forward, 
backward and central differences repeatedly. 

Theorem 4.1. Let T be nonuniform. Then 

P(Qf(tk)) ¥= Q(Pf(h)) for any applicable t k G T, 

where P and Q are any two different difference operations denoted by D + , D- or 
5. 

Proof. We only need to show the cases involving D + (D_f(t k )) and D + (Sf(t k )). 
According to (3.1) and (3.2), we have 



n-k 



m+i) - f(t k ) f(t k ) - /(t fc _x) 



h k h k -\ 



/hk 



h k -if{t k+ i) - (h k + h k -i)f(t k ) + u n 
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By the same token, 

7(*fc+i) - /(**) - /(tfc-i) 



L >-{ D +j{ t k)) = 7 



fc-1 



h k -!f(t k+1 ) - (h k + /ifc_i)/(t fc ) + h k f{t k _{) 



h kh\-i 

Recalling that h k / h k -i, we have 

D + (D_f(t k ))^D_(D + f(t k )). 

Further, 

'f(t k+2 )-f(t k ) f(t k+1 ) - /(i fe _!) 



//lfc-1 

(4.2) 



n c/rm ^ Wfc+i) ~ Wfc) 
D+{df{t k )) = — 



(h k + /»fc_i)(/(t fc +2) - /(**)) - (/ifc+i + h k )(f(t k +1 ) - /fa-i)) (4 3) 



5{D + f(t k )) = 



(h k+1 + h k )h k (h k + 
D+f(t k+1 ) - D+/(*fc-i) 



/(t fc+2 ) - /(ife+i) /(**) - /(tfc-i) 



/(/i fc + V 

= h k _i(f(t k+2 ) - f(t k+1 )) - h k+1 (f(t k ) - /fa-i)) 
h k +i(h k -i + h k )h k _i 

Therefore 



(4.4) 



S(D + f(t k ))jkD + (Sf(t k )) 
unless %_i = % = □ 
Corollary 4.1. Let T 6e uniform. Then we have 

P{Qf{tk)) = Q(Pf(t k )) for all applicable t k G T, 

where P and Q are any two different difference operations denoted by D + , D- or 
5. 

The proof of the corollary follows from (4.1)-(4.4). To explore interesting fea- 
tures of the proposed formulas on different meshes, we require exceptionally high 
orders of the derivative functions in the following theorem. The requirements can 
be conveniently eased when the types of the meshes are fixed. 

Theorem 4.2. Let T be nonuniform. Then none of the second order differences 
P(Qf(t k )), where P and Q are any two difference operations denoted by D + , D- 
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or 6, is a consistent approximation of f"(t k ). Further, if f is sufficiently smooth, 
D + (D + f(t k )) 



D_{D_f(t k )) 



S(6f(t k )) 



S(D+f(t k )) 



hk+i + h k „, . {hk+i + h k )(h k+ i + 2h k ) .,„. . 

"J W H tt; / (tk) 



+ 



2h k 
hk+i + hk 
24h k h k+1 



Qh k 

(% + i + %) 3 / (4) (l)-^/ (4) (6 



%-i + hk-2 j/, _ (hk-i + h k -2){2hk-i + hk-2) jw^ t ^ 



+ 



2hk~\ 
hk-i + 
24hk-ihk- 



2 L 



-^_i/ (4) (C) + (^-i + ^-2) 3 / (4) (C) 



(4.5) 



(4.6) 



2(h k + fc fc _x) ; W 

| (fefc+i + /tfc) 2 - (fefc_i + /lfc-2) 2 

6(h k + h k -i) 
(h k + /ife+i) 3 + {h k -i + /ifc-2) 3 



24(h k + h k -i) 
(h k+ i + M 4 / (5) (f) - (%_! + %_ 2 ) 4 / {5) (C) 



D + (D.f(t k )) 



£>_(£>+/(**)) 



D + (Sf(t k )) = 



120(h k + h fc _i) 
% + %-i ,///. % . h k — h\_ x 



+ 



2h k 
1 

120% 



-f"(t k ) + JL—hzlf" {tk ) + 
hif^(0-ht-J {5) (0 



120/i 



fc-i L 



^ + ^-1 
24% 

hi + ^-1 
6%_i J ' 24%_i 

^/ (5) (0-^-i/ (5) (C)l, 



hk + hk-i t iu \ h k - /i fc _ 1 
2fefc - J W H — / W + 

+ 



/ (4) (t fc ) 



f (A) (t k ) 



(4.7) 



(4.8) 



(4.9) 



%+i + h k -i „. . (%+i + %) 2 - /if. + frfchfc-i - /iLi ,////. >, 
"/ W H ttt / {tk) 



2h k 



+ 



1 



2\h k (h k + %_i) 

-^/ (4) (0-^-i/ (4) (0 



6% 

(% + %_!)(% +1 + %) 3 / {4) (|) 



(4.10) 



%+i + 2h k + %_i „ (h k + h k+ if ~h\- hk+ih 2 ^ 

—J V-k) H :rr 77 — ; r n J {tk) 



2{h k + h k . 



Qh k+ i(h k + h k . 



2Ah k+l (h k + h k -i) 
hUh k+1 f^(() 



(% + i + %)¥ 4) (0-^/ (4) (0 



(4.11) 
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D - mtk)) = ~^r ntk) + e^+^o^-i f {tk) 



i 



^/ (4) (0 + (^-i + ^-2)V (4) (C) 



24(h k + %_!)%_! 
+(/i fc + + %_ 2 ) 3 / (4) (C)] , (4.12) 

= 2(/> fc + /*_ 1 ) /( ^ ) + Mi^j f {tk) 

1 k- 2 ^/ (4) (o-^-i/ (4) (c) 



24(/i fc + h k -i)h k -2 
+(^-i + /» fc - 2 ) 4 / (4) (C)], (4-13) 

where £, £, £ and £ o^e different numbers and t k+2 > £ > t k , t k+ \ > £ > > 

C > *fc-l, t k > ( > t k -2- 

Proof. Due to their similarities, here we only present proofs of a few cases as an 
illustration. Let us consider D + (D+f(t k )), S(Sf(t k )), and D-(5f(t k )). It will be 
sufficient to show (4.5), (4.7) and (4.12). We have the Maclaurin series 

m+2) = f(t k ) + (h k+1 + h k )f\t k ) + (/tfc+i 2 + %) V fa) 

+**)y (ffc) + . . . + ( 4.i4) 

b n! 

/(tfc-2) = /fa) ~ + h k - 2 )f'(t k ) + f"(t k ) 

- (/lfc - 1+ 6 /lfc ' 2)3 /"' (t fc ) + • • • + (-ijn ^-i + ^-a)" /(»)(£). (4.15) 
Similar to (4.1) and (4.2), we have 



h k f{t k+ 2) - (h k+ i + h k )f( t k+ i) + h k+ if(t k ) 
h k +ihl 



D + (D + f(t k )) = \ Uk + 2 > _ v°k+i -r ^fe^/ yt-fc+i; t "-fc+u y-fc; ^ ^ 



Now, according to (3.7), (4.14), 
h k f{t k+ 2) - (h k+ i + h k )f(t k+ i) + h k+ if(t k ) 

= h k (f(t k) + + M/'fa) + {hk+1 + hk? f"{t k) + ...+ ( ^ +1 + fefc)n /W(o) 

- (h k+1 + %) (/(t fc ) + /*/'(<*) + -^/"(tfc) + • • • + ^/ (n) (0j + 



h k h k+1 (h k+1 + h k ) „ h k (h k+l + h k ) [{h k+1 + h k ) 2 - h\] 

-J ytk) H ~, / (tk) H 



3! 



+ 



h k (h k+ i + h k ) 



(%+i+%r- i / (n) (i)-^r 1 / (n) (o 



h k h k+ i(h k+ i + hk) , h k h k+1 (h k+ i + h k )(h k+ i + 2h k ) w 



3! 



/'"(**) + 



+ 



h k (h k+ i + /i A 



n! 



(/» fc+ i +/i fc ) n " 1 / (n) (i)-/»r 1 / (n) (o' 



Substituting the above into (4.16) and letting n = 4, we obtain (4.5) which indicates 
that the difference is not a consistent approximation of f"{t k ) unless h k+ \ = h k . Our 
next case involving (4.7) is interesting. To see this, we have 



W(t*)) 

7(4+2 



% + /i fc _i 

/(**) /(**) - f(t k - 2 ) 



/(h k + %_i) 



^fe+i + h k hk-i + /ifc_2 

_ (hk-i + h k - 2 )f(t k+2 ) - (hk+i + h k + hk-i + h k - 2 )f(t k ) + (fe fc +i + h k )f(t k - 2 ) 
(h k +i + /ifc)(/ifc + h k -i)(h k -i + %-2) 

We obtain (4.7) readily from the above by substituting (4.14) and (4.15). Further, 
using steps similar to those in the above discussion, we have 

D-(6f(t k )) = [(%-i + h k - 2 )f(t k+1 ) - (h k + h k -i)f(t k ) - (%_i + fc fc - 2 )/(tfc-i) 
+(h k + h k -i)f(t k - 2 )] /[(h k + h k -i)hk-i{h k -i + h k - 2 )]. 

Substitute (3.7), (3.8) and (4.15) into the above equation. We acquire (4.12) after 
simplification. Obviously, (4.12) indicates again that the finite difference is not a 
consistent approximation of f"(t k ) unless h k = h k -i = h k - 2 - □ 

Corollary 4.2. Let T be any mesh. Then 

£>+(£>_/(t fc )) = ^D.(D + f(t k )), t k € T. 
hk 



Corollary 4.3. Let a be a positive constant, and let ht+\ = aha for all possi- 
ble £ on a nonuniform mesh T. Then a necessary and sufficient condition for any 
aforementioned finite difference formula to be a consistent approximation of f" is 

a = 1. 
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Proof. Because the proofs of different cases are similar, we will show only one of 
them. Consider 5(D_f(tk)). To have a consistent approximation, we must require 
the coefficient of f"[tk) i n (4-13) to be one, that is, 

h k + Zhk-i + hk-2 _ a 2 h k _ 2 + 2ahk-2 + hk-2 _ a 2 + 2a + 1 _ 
2{h k + h k -i) ~ 2(a 2 h k ^ 2 + ah k _ 2 ) ~ 2(a 2 + a) 

Solving the above equation, we acquire a 2 + 2a + 1 = 2a 2 + 2a which implies that 
a 2 = 1. Since a is positive, therefore a = 1 is our only solution. □ 

Theorem 4.3. Let T be uniform. Then for any valid index k, 5(5 f(t k )), D + (D_f(t k )) 
and D_(D + f(t k )) are second order approximations of f"(t k ) if f is four times con- 
tinuously differentiable, and all other second order differences are first order approx- 
imations of f"(t k ) if f is three times continuously differentiable. Further, 

S(6f(t k ) 

D + (D_f(t k ) 
D-(D + f(t k ) 
D+(D + f{t k ) 

D+{5f(t k ) 
S(D+f(t k ) 
D-(6f(t k ) 
6{D-f(t k ) 

where £, £, ( and ( are constants for which t k+ 2 > £ > t k , t k+ \ > £ > t k , t k > ( > 
tk-i, tk > C > tk-2- 

Proof. We only need to show the nine equations listed. This can be implemented 
by letting h k+ i = h k = h k _\ = h k _ 2 = h in (4.5)-(4.13) and choosing proper values 
of n in the Maclaurin series (3.7), (3.8), (4.14) and (4.15). □ 

Remark 4.1. The difference formula (1.2) is a consistent first order approximation 
of f"(t) on a uniform mesh. 

Remark 4.2. The right-hand sides of the nine equations in Theorem 4.3 can be 
viewed as errors of the respective approximations on uniform meshes. 



/"(**) 

f"{t k ) 
f"(t k ) 
f"(t k ) 

f"(t k ) 

f"(t k ) 
f"{t k ) 
f"(t k ) 
f"(t k ) 



h 2 



/ (4) (l) + / (4) (0 

/ (4) (0+/ (4) (C) 



6 

lf_ r 

24 

D + (D.f(t k ))-f"(t k ), 
h 



4/ w (o - no 
no-4no 
sno-no 



no 



3 
h 

3 
h 

12 

D + (5f(t k ))-f"(t k ), 

y 2 [no + /"'(c) - sno 

D-(5f(t k ))-f"(t k ), 
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Remark 4.3. Although none of the standard second order differences discussed in 
this paper may approximate / on a nonuniform mesh, specially formulated finite 
difference formulas may work. An example is 

(h k -.i + h k )/2 

It can be shown by using Maclaurin series that (4.17) provides a first order approx- 
imation of f"(tk) on any nonuniform mesh. 



5. Numerical examples. It has been known that trigonometric functions 
y = sin(crf), y = cos(bt) possess excellent smoothness for applications. The functions 
are infinitely differentiable. To illustrate our results, we choose the function 



y = -sin(47rt), < t < 1. 
The corresponding derivatives are 



y 



-47rcos(47rt), y" = (4vr) 2 sin(4vrt), < t < 1. 



(5.1) 



(5.2) 



uniform & nonuniform meshes 



uniform & nonuniform meshes 




Figure 5.1. Plot of grid step sizes (red dotted line is for the uniform mesh and blue solid 
curve and area are for the nonuniform mesh, LEFT) . In this figure, [0, 1] is divided into 22 
subregions proportional to the 22 steps used. Height of each of the bars located in the subregion 
is the actual size of the corresponding step size; The second figure is for smoothness ratio for 
the uniform mesh (red dotted line) and the nonuniform mesh (blue curve). In the case, [0, 1] is 
divided into 21 subregions proportional to the first 21 steps used. Height of each bar is for the 
corresponding ratio. 

To achieve better simulation resolution, we add eleven additional external points 
into the uniform and nonuniform meshes introduced in Example 2.1, respectively. 
The new points are separated by the original points. The new points in the nonuni- 
form mesh are particularly set to be closer to their right-side points, as plot- 
ted in Figure 5.1 (LEFT). In the second frame of Figure 5.1, we show the ratio 
hk+i/hk, k = 1, 2, . . . , 21, (the red line is for the uniform mesh and the blue curve 



12 



is for the nonuniform mesh). A logarithmic scale in the Y-direction is used to show 
more precisely the ratio, which is often called the smoothness ratio in engineering 
computations and is chosen between and 10. 

Numerical computations will be carried out for a number of the differences. In 
our experiments, a numerical solution g is acceptable if sgei(/, g) <fc<Cl, and it is 
unacceptable if sgei(/, g) > 1. Matlab® will be used throughout the experiments. 

Example 5.1. Consider the first difference 5f. The true first derivative in (5.2), 
and the central difference are plotted in Figure 5.2. Though the approximations in 
both cases are consistent, we may observe that the approximation in the uniform 
case is much better than that in the nonuniform case, since in the former situation 
the approximation is second order. 



uniform mesh case nonuniform mesh case green dash-uniform: black-nonuniform 




t t t 



Figure 5.2. Plot of true solution (red dotted curves) and approximations (blue curves) on the 
uniform mesh (LEFT) and the nonuniform mesh (MIDDLE). Corresponding scaled local differ- 
ence curves are given (RIGHT: green for the uniform mesh case and black for the nonuniform 
case). We may observe from the first two pictures that, even though approximations on both 
meshes are consistent, the approximation on the uniform mesh is much "nicer" than that on the 
nonuniform mesh. The third picture confirms this by showing the scaled difference values over 
the domain (the black curve acts more violently with a relatively large amplitude). 

These are further confirmed by the third frame in Figure 5.2 in which relative 
errors are given. While the scaled global error indicator in the uniform mesh case 
sgei ~ 0.0535, in the nonuniform case, it reaches 0.4759 which is about 9 times 
the uniform mesh case! The irregular smoothness ratios shown in Figure 5.1 may 
explain why the error in the nonuniform grid case is oscillatory. 

Example 5.2. Consider the second difference D + (D + f). We plot the true sec- 
ond derivative given in (5.2) and the difference in Figure 5.3. It can be observed 
that while the difference approximates the derivative function reasonably over the 
uniform mesh (global error indicator value sgei ~ 0.5559 because relatively large 
steps are used), it produces an unacceptable results on the nonuniform mesh with 
an indicator value sgei ~ 2.5487. The irregular oscillations of the finite difference on 
nonuniform mesh is also unacceptable. 
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uniform mesh case nonuniform mesh case green dash-uniform; black-nonuniform 




t 06 ° 7 ° 8 ° 9 ° 1 ° 3 ° 3 34 t 36 ° 7 °' 3 33 ° 31 3 ' 3 33 °' 4 t 



Figure 5.3. Plot of true solution (red dotted curves) and approximations (blue curves) on the 
uniform mesh (LEFT) and the nonuniform mesh (MIDDLE). Corresponding scaled local differ- 
ence curves are given (RIGHT: green for the uniform mesh case and black for the nonuniform 
case). We may observe from the first two pictures that, while the approximation on the uniform 
mesh is still reasonable, the approximation on the nonuniform mesh becomes unacceptable. The 
third picture confirms these by showing the scaled difference values over the domain (the black 
curve acts violently with a significantly large amplitude). 

Example 5.3. Consider the second difference 5{D_f). The intention of using the 
central difference is to improve the numerical result. It is found in Figure 5.4 that 

sgei ~ 0.2817 when T is uniform 
sgei ~ 0.7829 when T is nonuniform 

In addition to a larger global relative error, in Figure 5.4, we may also observe that 
if a nonuniform mesh is used, there is a great irregularity in errors. If the difference 
is used as an approximation, very incorrect answers will probably result. 



uniform mesh case nonuniform mesh case green dash-uniform; black-nonuniform 




t t t 



Figure 5.4. Plot of true solution (red dotted curves) and approximations (blue curves) on the 
uniform mesh (LEFT) and the nonuniform mesh (MIDDLE). Corresponding scaled local differ- 
ence curves are given (RIGHT: green for the uniform mesh case and black for the nonuniform 
case). We may observe from the first two pictures that, while the approximation on the uniform 
mesh looks good, the approximation on the nonuniform mesh becomes unacceptable. The third 
picture again confirms these by showing the scaled difference values over the domain (the black 
curve oscillates violently with a large amplitude). 
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Example 5.4. Consider the second difference 5(5 f). Although the use of central 
differences continues to improve the numerical result, it still cannot change the basic 
features of the approximations. In this case, we have 

sgei ~ 0.1031 when T is uniform 
sgei ~ 0.5080 when T is nonuniform 

The strong irregularity in the nonuniform mesh case shown in Figure 5.5 is not 
surprising. 

uniform mesh case nonuniform mesh case green dash-uniform; black-nonuniform 



Figure 5.5. Plot of true solution (red dotted curves) and approximations (blue curves) on the 
uniform mesh (LEFT) and the nonuniform mesh (MIDDLE). Corresponding scaled local differ- 
ence curves are given (RIGHT: green for the uniform mesh case and black for the nonuniform 
case). We may observe from the first two pictures that, while the approximation on the uniform 
mesh is improved due to the use of central difference, the inconsistency of the formula on the 
nonuniform mesh remains the same. These are again confirmed by the third picture through 
the scaled difference values over the domain (the black curve oscillates successively with a large 
amplitude) . 



Example 5.5. Let us consider the simple harmonic oscillator problem where neither 
a driving force nor friction is assumed [9]. If 4>(t) is the displacement of the system 
at time t, then the second derivative 4>"(t) is its acceleration. Based on Hooke's 
Law and Newton's Second Law, we obtain the following second order differential 
equation, 

<t>"(t) = -K<t>(t), t > t , (5.3) 
together with the initial conditions 

<K*o) = l, 0'(*o) = -l, (5.4) 

where k is a positive constant. 

It is not difficult to verify [10] that the solution of (5.3), (5.4) is 

4>(t) = — — sin yfnt + cos \fnt, t>to. (5.5) 
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On the other hand, replacing the second derivative in (5.3) by the backward-forward 
difference, we acquire 

D_{D + w) = -Kw(t), teT, (5.6) 

where T is a mesh over the interval [to, b]. We wish the solution of (5.6), (5.4) to be 
an approximation of the solution of (5.3), (5.4). 

First, we let Ti be a nonuniform mesh over the interval [to,b] with decreasing 
steps /ifc = rhk-i, k = 1, 2, . . . , m, for ho = 1/10, r = 50/59 and m = 200. Second, 
we let T2 be an uniform mesh with h = (b — to)/10. 

m 

Let k = 4vr 2 , t = and b =J2 h k ~ 59/90. In Figure 5.6 (LEFT), we plot the 

fc=0 

numerical solutions of (5.6), (5.4) on Ti (blue dashed curve), exact solution of the 
initial value problem (5.3), (5.4) (red dotted curve) together with a numerical solu- 
tion of (5.6), (5.4) on T2 (green curve). A forward difference is used to approximate 
the first derivative in (5.4). Although errors defined for the numerical solution of 
initial value problems are slightly different from those used for function approxima- 
tions, to illustrate the inconsistency of second order derivative approximations on 
nonuniform meshes, let us continue to use the measurements introduced by Defi- 
nition 2.4. Scaled local differences of the numerical solutions on Ti (blue dashed 
curve) and on T2 (green curve) are given in Figure 5.6 (RIGHT). 



numerical and true solutions adjusted local differences of the solutions 




Figure 5.6. Plots of the solutions of the difference equation problem (5.6), (5.4) and the true 
solution (LEFT); and the scaled local difference values of the numerical solutions (RIGHT) on 
the nonuniform mesh and uniform mesh, respectively. While the numerical solution and related 
sld values on the nonuniform mesh Ti are given by blue dashed curves, the numerical solution 
and related sld values on the uniform mesh T2 are represented by solid green curves. We may 
observe in the first picture that, while the numerical solution on the uniform mesh maintains 
a good match to the true solution, which is indicated by the red dotted curve, the numerical 
solution on the nonuniform mesh swings away as t increases. The phenomena are confirmed by 
the second picture through the scaled difference values over the domain used (the blue curve 
oscillates and its amplitude increases rapidly as t — >■ b). 

It is interesting to observe that while the numerical solution obtained on the 
uniform mesh keeps a steady error level from the true solution, although relatively 
large step is being used (h = 59/900 0.0656), the numerical solution obtained on 
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the nonuniform mesh runs away rapidly from the true solution, despite of the fact 
that finer steps are employed (hk < 0.001 for k > 14). 

The global error indicator values for the two cases are also significantly different. 
It is found that sgei « 0.5055 on Ti and occurs at x = b while sgei « 0.1945 on T2 
and occurs at x ~ 0.2622. 

Since the sizes of h and do not play a major role in the aforementioned error 
phenomena, what can be the main cause for the unsatisfactory approximation on 
Ti? 

To see the answer, we may notice that 

hk + hk-i = r + l ^ 1 
2hk-i 2 

Therefore, according to (4.9) in Theorem 4.2, difference equation (5.6) actually 
approximates 

-^V(t) = -K^(t), t > t , 

instead of (5.3)! Thus the numerical solution cannot be satisfactory. Of course, 
there are several factors, such as stability, that may affect the numerical solution. 
Needless to say, however, consistency is the most fundamental factor. 

Conversely, it is not hard to verify that the improved T>2 difference given in 
(4.17) offers a good approximation. The reader may wish to experiment with the 
interesting computations! 



6. Conclusions. From the foregoing discussions and numerical experiments, 
we may conclude that: 

1. The key issue for any finite difference formula is its consistency. Only con- 
sistent formulas can be used for approximating derivatives. The effectiveness 
of an approximation can be measured by its order of accuracy. For instance, 
(1.1) is an order one formula according to Theorems 3.1 and 3.2. 

2. Consistent finite difference approximations can be derived on uniform or nonuni- 
form grids. Generally speaking, the higher the order of accuracy is, the more 
accurate the finite difference formula can be. 

3. Difference approximations on uniform meshes cannot be applied blindly on 
nonuniform meshes, nor can difference formulas be composed to form consis- 
tent approximations to second derivatives. At best, they may lose accuracy; 
at worst they are not consistent. 

4. Assuming that the functions involved are sufficiently smooth, errors of ap- 
proximations can be estimated by using Maclaurin series. Errors in numerical 
experiments can also be computed via local and global error formulas. 
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5. The consistency and accuracy of different approximations can be demonstrated 
through the use of computer simulations. Since simulations are based on 
particular examples, they are not as rigorous as mathematical proofs. 
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